

clear all
clc

% b = str2double(getenv('SLURM_ARRAY_TASK_ID'));

b       = 2;
v       = 2;
datsize = 100;
S       = 200; B = 100;
out     = 1;

load(['DataResult_v' num2str(v) '_ival' num2str(b) '_myopic_out_' num2str(datsize) '.mat'])

theta = theta_score;

%% Halton draws (use a different set of S draws for each MC)

StateVec = randi(1e9,[B,2,S]);

FirmError_Sim = zeros(N,S);
ParcelError_Sim = zeros(N,S);

for s = 1:S

const.NObs = N; 

rng(StateVec(v,1,s))
skip = randi(1e3);
leap = randi(1e2);

p = haltonset(1,'Skip',skip,'Leap',leap);
p = scramble(p,'RR2');
draw = net(p,const.NObs);
draw = norminv(draw);

% FirmError_Sim(:,s) = draw(firmID,:);
FirmError_Sim(:,s) = draw;

clear draw p

const.NObs = N; 

rng(StateVec(v,2,s))
skip = randi(1e3);
leap = randi(1e2);

p = haltonset(1,'Skip',skip,'Leap',leap);
p = scramble(p,'RR2');
draw = net(p,const.NObs);
draw = norminv(draw);

ParcelError_Sim(:,s) = draw;

clear draw p

end

%% Verify Equilibrium Assignment

% Parcel preferences
gamma = [1 theta(1,1:NN - 1)];
% Firm preferences
delta = [1 theta(1,NN:NN + JJ - 2)];
alpha = theta(1,NN + JJ - 1);

% Firm Value of parcels
ParcelUtility = ParcelIndex*gamma';
% Parcel value of firms (and leases)
FirmUtility = FirmMatrix*delta';

shareOLD_st = share0;
iter = 0;
tol = 10;
while tol > 1e-6 && iter < 300

    % Sub-market profits and index
    Offer_sub  = ParcelUtility + alpha*shareOLD_st; 
    Accept_sub = FirmUtility; 

    Aoffer  = reshape(Offer_sub,[J,N/J]); 
    Baccept = reshape(Accept_sub,[J,N/J]);

    % Sort preferences
    if out == 1 
        rankA = SortingMatInC_out(Aoffer); 
        rankB = SortingResMatInC_out(Baccept);
    else
        rankA = SortingMatInC(Aoffer); 
        rankB = SortingResMatInC(Baccept);
    end
    [out,FirmMatchOUT]      = MatchAofferBaccept_GS_out(rankA,rankB,FirmMatchConstraint);
    [~,shareNEW_st]  = ShareMarket_st(FirmMatchOUT',MarketCount0,J);

    tol = sum(abs(shareOLD_st - shareNEW_st));
    shareOLD_st = shareNEW_st;

%    shareStOUT_org = (1/2)*(shareStOUT + shareStOUT_org);
    iter = iter + 1;
end

%% Prepare Data Moments for numerical derivative

% StateVec = randi(1e9,[NN + JJ + 5,3,S]);
% rng(StateVec(b,1,s))

% s_rand = StateVec(1,3,1);
% s_rand = RandStream('mcg16807','Seed',s_rand);
% RandStream.setDefaultStream(s_rand);

% FirmError_Sim = randn(N,S);
% ParcelError_Sim = randn(N,S);

match.ParcelError_Sim   = ParcelError_Sim;
match.FirmError_Sim     = FirmError_Sim;

match.MarketIndex       = MarketIndex;

match.MarketCount       = MarketCount;
match.FirmMatch         = FirmMatch;
match.ParcelIndex       = ParcelIndex;
match.FirmMatrix        = FirmMatrix;

match.MarketIndex0      = MarketIndex0;
match.FirmMatch0        = FirmMatch0;
match.ParcelIndex0      = ParcelIndex0;
match.FirmMatrix0       = FirmMatrix0;

match.NN                = NN;
match.JJ                = JJ;
match.N                 = N;
match.J                 = J;
match.M                 = M;
match.S                 = S;

match.out               = out;

match.FirmShare0        = FirmShare0;             % Share by market (J*M x 1)
match.share0            = share0;                 % Stacked share (J*N x 1)

match.FirmMatchConstraint   = FirmMatchConstraint;
match.MarketCount0          = MarketCount0;

MomentSwitch = [1 1 1 1 1 ];
match.MomentSwitch = MomentSwitch;

[PMean0,FMean0,JntDist0,Cov0,Pvar0,Fvar0,PFcovar0,Pmean_firm0,Pvar_firm0,PFJntDist_firm0,PFcovar_firm0,~,~,~,~] = GroupMoments0_new(match.ParcelIndex0,match.FirmMatrix0,match.FirmMatch0,J,MarketIndex0,MarketCount0,MomentSwitch);
[GpPMeanShare0,GpFMeanShare0,GpShareMarket0,GpShareFirm0,GpFMean0,GpParcelCov0,GpFirmCov0] = GroupMoments0_share2(match.ParcelIndex0,match.FirmMatrix0,match.FirmMatch0,J,MarketIndex0,MarketCount0,FirmShare0);

temp                    = [FMean0 PMean0];
match.MMean_Moment      = temp(:);
match.MCov_Moment       = PFcovar0(:);
temp                    = [Fvar0 Pvar0];
match.MVar_Moment       = temp(:);

match.JntDist0          = JntDist0(:);
match.Cov0              = Cov0(:);

match.JMean_Moment      = Pmean_firm0(:);
match.JCov_Moment       = PFcovar_firm0(:);
match.JVar_Moment       = Pvar_firm0(:);
match.JDist_Moment      = PFJntDist_firm0(:);

% match.meanGroup_Moment  = GpMean0(:);
% match.varGroup_Moment   = GpVar0(:);

% match.JMVar_Moment      = GpVar_byJ0(:);
% match.JMMean_Moment     = GpMean_byJ0(:);

temp                    = [GpFMeanShare0 GpPMeanShare0];
match.GpMeanShare0      = temp(:);
match.GpShareMarket0    = GpShareMarket0;
match.GpShareFirm0      = GpShareFirm0;
match.ShareFMean0       = GpFMean0;
temp                    = [GpFirmCov0 GpParcelCov0];
match.GpCovShare0       = temp(:);

%% J-stat

[M_sim,M_data,J_stat,Share_Sim,critical_value] = Jstat_end_myopic(betaStar,match);

disp 'J-stat'
J_stat
disp 'cannot reject the null that the moments fit if J is less than crit'
disp 'critical val'
critical_value

%% Data for the observed match

h = 1e-2;

betaStar_m = theta;
betaStar_p = theta;

M_sim_m = zeros(1,JJ + NN); M_m = zeros(1,JJ + NN);
M_sim_p = zeros(1,JJ + NN); M_p = zeros(1,JJ + NN);

for b = 1:size(betaStar,2)
    betaStar_m(1,b) = betaStar(1,b) - h;

    % Parcel preferences
    gamma = [1 betaStar_m(1,1:NN - 1)];
    % Firm preferences
    delta = [1 betaStar_m(1,NN:NN + JJ - 2)];
    alpha = betaStar_m(1,NN + JJ - 1);

    % Firm Value of parcels
    ParcelUtility = ParcelIndex*gamma';
    % Parcel value of firms (and leases)
    FirmUtility = FirmMatrix*delta';
    
    SampleMMean_SimMoment       = zeros(M,JJ + NN,S);
    SampleMVar_SimMoment        = zeros(M,JJ + NN,S);
    SampleMCov_SimMoment        = zeros(M,JJ*NN,S);

    SampleJMean_SimMoment       = zeros(J,NN,S);
    SampleJCov_SimMoment        = zeros(J,NN*JJ,S);
    SampleJVar_SimMoment        = zeros(J,NN,S);
    SampleJDist_SimMoment       = zeros(J,NN*JJ,S);

    SampleJntDist_SimMoment     = zeros(NN*JJ,S);
    SampleCov_SimMoment         = zeros(NN*JJ,S);

    Share_Sim                   = zeros(M*J,S);
    ShareMean_Sim               = zeros(M,NN + JJ,S);
    ShareFirm_Sim               = zeros(M,S);
    ShareCov_Sim                = zeros(M,NN + JJ,S);

    parfor s = 1:S
%        rng(StateVec(b,1,s))
%        s_rand = StateVec(b,1,s);
%        s_rand = RandStream('mcg16807','Seed',s_rand);
%        RandStream.setDefaultStream(s_rand);

%        FirmError_Sim = randn(N,1);
%        ParcelError_Sim = randn(N,1);

%        FirmMatchFinal = zeros(1,1);
%        shareOUTFinal = zeros(1,1);

        shareIter = share0;
        shareKEEP = []; shareStKEEP = [];
        tol = 1; iter = 1;
        while tol > 1e-6 && iter < 300
            % Sub-market profits and index
            ParcelUtility_sub = ParcelUtility + alpha*shareIter + ParcelError_Sim(:,s); 
            FirmUtility_sub = FirmUtility + FirmError_Sim(:,s); 

            Aoffer = reshape(ParcelUtility_sub,[J,N/J]); 
            Baccept = reshape(FirmUtility_sub,[J,N/J]);

            % Sort preferences
            if out == 1
                rankA = SortingMatInC_out(Aoffer); 
                rankB = SortingResMatInC_out(Baccept);
            else
                rankA = SortingMatInC(Aoffer); 
                rankB = SortingResMatInC(Baccept);
            end
            
            % GS algorithm
            [~,FirmMatchOUT] = MatchAofferBaccept_GS_out(rankA,rankB,FirmMatchConstraint);
            [shareOUT,shareStOUT] = ShareMarket_st(FirmMatchOUT',MarketCount0,J);
            
            shareStKEEP = [shareStKEEP shareStOUT];
            shareKEEP   = [shareKEEP shareOUT];
            
            % tol = mean(abs(share0 - shareStOUT));
            if iter > 1
            %   tol = mean(abs(shareKEEP(:,iter - 1) - shareKEEP(:,iter)))
                tol = sum(abs(shareKEEP(:,iter - 1) - shareKEEP(:,iter)));
            else
            %   tol = mean(abs(share0 - shareStOUT))
                tol = sum(abs(shareIter - shareStOUT));
            end
            shareIter = shareStOUT;
%            shareIter = mean(shareStKEEP,2);
%            FirmMatchFinal(1,1:size(FirmMatchOUT,2)) = FirmMatchOUT;
%            shareOUTFinal(1:size(shareOUT,1),1) = shareOUT;
            iter = iter + 1;
        end
        
        % Firm group royalty/clause means for moments
        [~,~,GpPMeanShare,GpFMeanShare,~,GpShareFirm,~,GpParcelCov,GpFirmCov] = GroupMoments_share2(ParcelIndex,FirmMatrix,FirmMatchOUT,J,MarketIndex0,MarketCount0,shareOUT);
        [~,~,PMean,FMean,JntDist,Cov,Pvar,Fvar,PFcovar,Pmean_firm,Pvar_firm,PFJntDist_firm,PFcovar_firm,~,~,~,~] = GroupMomentsSmall_new(ParcelIndex,FirmMatrix,FirmMatchOUT,J,MarketIndex0,MarketCount0,MomentSwitch);

        % Calculate the moments
        SampleMMean_SimMoment(:,:,s) = [FMean PMean];
        SampleMCov_SimMoment(:,:,s)  = PFcovar;
        SampleMVar_SimMoment(:,:,s)  = [Fvar Pvar];

        SampleJMean_SimMoment(:,:,s) = Pmean_firm;
        SampleJCov_SimMoment(:,:,s)  = PFcovar_firm;
        SampleJVar_SimMoment(:,:,s)  = Pvar_firm;
        SampleJDist_SimMoment(:,:,s) = PFJntDist_firm;

        SampleJntDist_SimMoment(:,s) = JntDist;
        SampleCov_SimMoment(:,s)     = Cov;

        % Calculate share - Market
        Share_Sim(:,s)          = shareOUT;
        ShareMean_Sim(:,:,s)    = [GpFMeanShare GpPMeanShare];
        ShareFirm_Sim(:,s)      = GpShareFirm;
        ShareCov_Sim(:,:,s)     = [GpFirmCov GpParcelCov];
    
    end
    
    JntDist_SimMoment    = sum(SampleJntDist_SimMoment,2)*(1/S);
    Cov_SimMoment        = sum(SampleCov_SimMoment,2)*(1/S);

    % Market Moments
    MMean_SimMoment      = sum(SampleMMean_SimMoment,3)*(1/S);
    MCov_SimMoment       = sum(SampleMCov_SimMoment,3)*(1/S);
    MVar_SimMoment       = sum(SampleMVar_SimMoment,3)*(1/S);

    % Firm Moments
    JMean_SimMoment      = sum(SampleJMean_SimMoment,3)*(1/S);
    JCov_SimMoment       = sum(SampleJCov_SimMoment,3)*(1/S);
    JVar_SimMoment       = sum(SampleJVar_SimMoment,3)*(1/S);
    JDist_SimMoment      = sum(SampleJDist_SimMoment,3)*(1/S);

    % Market Structure Moments
    Share_Sim           = sum(Share_Sim,2)*(1/S);
    ShareMean_Sim       = sum(ShareMean_Sim,3)*(1/S);
    ShareFirm_Sim       = sum(ShareFirm_Sim,2)*(1/S);
    ShareCov_Sim        = sum(ShareCov_Sim,3)*(1/S);
    
    M_sim = [JntDist_SimMoment; Cov_SimMoment; ...
        MMean_SimMoment(:); MCov_SimMoment(:); MVar_SimMoment(:); ...
        JMean_SimMoment(:); JCov_SimMoment(:); JVar_SimMoment(:); JDist_SimMoment(:)];
    M_data = [match.JntDist0; match.Cov0; ...
        match.MMean_Moment; match.MCov_Moment; match.MVar_Moment; ...
        match.JMean_Moment; match.JCov_Moment; match.JVar_Moment; match.JDist_Moment];

    % Full share set
    M_share_sim = [Share_Sim; ShareMean_Sim(:); ShareFirm_Sim; ShareCov_Sim(:)];
    M_share = [match.FirmShare0; match.GpMeanShare0; match.GpShareFirm0; match.GpCovShare0];

    M_sim_m(1:size(M_sim,1) + size(M_share_sim,1),b) = [M_sim; M_share_sim];
    M_m(1:size(M_sim,1) + size(M_share_sim,1),b) = [M_data; M_share];
end

for b = 1:size(betaStar,2)
    betaStar_p(1,b) = betaStar(1,b) + h;

    % Parcel preferences
    gamma = [1 betaStar_p(1,1:NN - 1)];
    % Firm preferences
    delta = [1 betaStar_p(1,NN:NN + JJ - 2)];
    alpha = betaStar_p(1,NN + JJ - 1);

    % Firm Value of parcels
    ParcelUtility = ParcelIndex*gamma';
    % Parcel value of firms (and leases)
    FirmUtility = FirmMatrix*delta';

    SampleMMean_SimMoment       = zeros(M,JJ + NN,S);
    SampleMVar_SimMoment        = zeros(M,JJ + NN,S);
    SampleMCov_SimMoment        = zeros(M,JJ*NN,S);

    SampleJMean_SimMoment       = zeros(J,NN,S);
    SampleJCov_SimMoment        = zeros(J,NN*JJ,S);
    SampleJVar_SimMoment        = zeros(J,NN,S);
    SampleJDist_SimMoment       = zeros(J,NN*JJ,S);

    SampleJntDist_SimMoment     = zeros(NN*JJ,S);
    SampleCov_SimMoment         = zeros(NN*JJ,S);

    Share_Sim                   = zeros(M*J,S);
    ShareMean_Sim               = zeros(M,NN + JJ,S);
    ShareFirm_Sim               = zeros(M,S);
    ShareCov_Sim                = zeros(M,NN + JJ,S);

    parfor s = 1:S
%         rng(StateVec(b,2,s))
%        s_rand = StateVec(b,2,s);
%        s_rand = RandStream('mcg16807','Seed',s_rand);
%        RandStream.setDefaultStream(s_rand);

%        FirmError_Sim = randn(N,1);
%        ParcelError_Sim = randn(N,1);

        shareIter = share0;
        shareKEEP = []; shareStKEEP = [];
        tol = 1; iter = 1;
        while tol > 1e-6 && iter < 300
            % Sub-market profits and index
            ParcelUtility_sub = ParcelUtility + alpha*shareIter + ParcelError_Sim(:,s); 
            FirmUtility_sub = FirmUtility + FirmError_Sim(:,s); 

            Aoffer = reshape(ParcelUtility_sub,[J,N/J]); 
            Baccept = reshape(FirmUtility_sub,[J,N/J]);

            % Sort preferences
            if out == 1
                rankA = SortingMatInC_out(Aoffer); 
                rankB = SortingResMatInC_out(Baccept);
            else
                rankA = SortingMatInC(Aoffer); 
                rankB = SortingResMatInC(Baccept);
            end
            
            % GS algorithm
            [~,FirmMatchOUT] = MatchAofferBaccept_GS_out(rankA,rankB,FirmMatchConstraint);
            [shareOUT,shareStOUT] = ShareMarket_st(FirmMatchOUT',MarketCount0,J);
            
            shareStKEEP = [shareStKEEP shareStOUT];
            shareKEEP   = [shareKEEP shareOUT];
            
            % tol = mean(abs(share0 - shareStOUT));
            if iter > 1
            %   tol = mean(abs(shareKEEP(:,iter - 1) - shareKEEP(:,iter)))
                tol = sum(abs(shareKEEP(:,iter - 1) - shareKEEP(:,iter)));
            else
            %   tol = mean(abs(share0 - shareStOUT))
                tol = sum(abs(shareIter - shareStOUT));
            end
            shareIter = shareStOUT;
%            shareIter = mean(shareStKEEP,2);
%             FirmMatchFinal(s,1:size(FirmMatchOUT,2)) = FirmMatchOUT;
%             shareOUTFinal(1:size(shareOUT,1),s) = shareOUT;

            iter = iter + 1;
        end
        
        % Firm group royalty/clause means for moments
        [~,~,GpPMeanShare,GpFMeanShare,~,GpShareFirm,~,GpParcelCov,GpFirmCov] = GroupMoments_share2(ParcelIndex,FirmMatrix,FirmMatchOUT,J,MarketIndex0,MarketCount0,shareOUT);
        [~,~,PMean,FMean,JntDist,Cov,Pvar,Fvar,PFcovar,Pmean_firm,Pvar_firm,PFJntDist_firm,PFcovar_firm,~,~,~,~] = GroupMomentsSmall_new(ParcelIndex,FirmMatrix,FirmMatchOUT,J,MarketIndex0,MarketCount0,MomentSwitch);

        % Calculate the moments
        SampleMMean_SimMoment(:,:,s) = [FMean PMean];
        SampleMCov_SimMoment(:,:,s)  = PFcovar;
        SampleMVar_SimMoment(:,:,s)  = [Fvar Pvar];

        SampleJMean_SimMoment(:,:,s) = Pmean_firm;
        SampleJCov_SimMoment(:,:,s)  = PFcovar_firm;
        SampleJVar_SimMoment(:,:,s)  = Pvar_firm;
        SampleJDist_SimMoment(:,:,s) = PFJntDist_firm;

        SampleJntDist_SimMoment(:,s) = JntDist;
        SampleCov_SimMoment(:,s)     = Cov;

        % Calculate share - Market
        Share_Sim(:,s)          = shareOUT;
        ShareMean_Sim(:,:,s)    = [GpFMeanShare GpPMeanShare];
        ShareFirm_Sim(:,s)      = GpShareFirm;
        ShareCov_Sim(:,:,s)     = [GpFirmCov GpParcelCov];    
    end
    
    JntDist_SimMoment    = sum(SampleJntDist_SimMoment,2)*(1/S);
    Cov_SimMoment        = sum(SampleCov_SimMoment,2)*(1/S);

    % Market Moments
    MMean_SimMoment      = sum(SampleMMean_SimMoment,3)*(1/S);
    MCov_SimMoment       = sum(SampleMCov_SimMoment,3)*(1/S);
    MVar_SimMoment       = sum(SampleMVar_SimMoment,3)*(1/S);

    % Firm Moments
    JMean_SimMoment      = sum(SampleJMean_SimMoment,3)*(1/S);
    JCov_SimMoment       = sum(SampleJCov_SimMoment,3)*(1/S);
    JVar_SimMoment       = sum(SampleJVar_SimMoment,3)*(1/S);
    JDist_SimMoment      = sum(SampleJDist_SimMoment,3)*(1/S);

    % Market Structure Moments
    Share_Sim           = sum(Share_Sim,2)*(1/S);
    ShareMean_Sim       = sum(ShareMean_Sim,3)*(1/S);
    ShareFirm_Sim       = sum(ShareFirm_Sim,2)*(1/S);
    ShareCov_Sim        = sum(ShareCov_Sim,3)*(1/S);
    
    M_sim = [JntDist_SimMoment; Cov_SimMoment; ...
        MMean_SimMoment(:); MCov_SimMoment(:); MVar_SimMoment(:); ...
        JMean_SimMoment(:); JCov_SimMoment(:); JVar_SimMoment(:); JDist_SimMoment(:)];
    M_data = [match.JntDist0; match.Cov0; ...
        match.MMean_Moment; match.MCov_Moment; match.MVar_Moment; ...
        match.JMean_Moment; match.JCov_Moment; match.JVar_Moment; match.JDist_Moment];

    % Full share set
    M_share_sim = [Share_Sim; ShareMean_Sim(:); ShareFirm_Sim; ShareCov_Sim(:)];
    M_share = [match.FirmShare0; match.GpMeanShare0; match.GpShareFirm0; match.GpCovShare0];

    M_sim_p(1:size(M_sim,1) + size(M_share_sim,1),b) = [M_sim; M_share_sim];
    M_p(1:size(M_sim,1) + size(M_share_sim,1),b) = [M_data; M_share];
end

%%

clear ParcelError_Sim FirmError_Sim
clear match

save(['DataResult_v' num2str(v) '_ival' num2str(b) '_share_' num2str(datsize) '_inf.mat'])

%%

% q moments x p parameters
Full_M_p = M_sim_p(:,1:NN+JJ-1);
Full_M_m = M_sim_m(:,1:NN+JJ-1);

M_share_prime = (Full_M_p - Full_M_m)./repmat(([betaStar] + h - ([betaStar] - h)),[size(Full_M_p,1),1]);
W = size(Full_M_p,1);
lhs_share = inv(M_share_prime'*W*M_share_prime)*M_share_prime'*W;
rhs_share = W'*M_share_prime*inv(M_share_prime'*W*M_share_prime);

Q_share_s = [[betaStar]' sqrt(diag(inv(M_share_prime'*W*M_share_prime)*(1 + 1/S)))]

Sigma_share = lhs_share*W*rhs_share;
StdError_share = sqrt(diag(inv(Sigma_share)));

bound = size(M_sim_p,1) - size(M_share,1);

M_prime = (M_sim_p(1:bound,1:NN+JJ-1) - M_sim_m(1:bound,1:NN+JJ-1))./repmat((betaStar + h - (betaStar - h)),[bound,1]);
W = eye(size(M_prime,1));
lhs = inv(M_prime'*W*M_prime)*M_prime'*W;
rhs = W'*M_prime*inv(M_prime'*W*M_prime);

Q_s = [betaStar' sqrt(diag(inv(M_prime'*W*M_prime)*(1 + 1/S)))]

Sigma = lhs*W*rhs;
StdError = sqrt(diag(inv(Sigma)));

save(['DataResult_v' num2str(v) '_ival' num2str(b) '_share_' num2str(datsize) '_inf.mat'])

quit
